home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qmomo.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-06-05  |  4.3 KB  |  194 lines
  1. /* integration/qmomo.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_integration.h>
  23. #include <gsl/gsl_errno.h>
  24.  
  25. static void
  26. initialise (double * ri, double * rj, double * rg, double * rh,
  27.         double alpha, double beta);
  28.  
  29. gsl_integration_qaws_table * 
  30. gsl_integration_qaws_table_alloc (double alpha, double beta, int mu, int nu)
  31. {
  32.   gsl_integration_qaws_table * t;
  33.  
  34.   if (alpha < -1.0)
  35.     {
  36.       GSL_ERROR_VAL ("alpha must be greater than -1.0", GSL_EINVAL, 0);
  37.     }
  38.  
  39.   if (beta < -1.0)
  40.     {
  41.       GSL_ERROR_VAL ("beta must be greater than -1.0", GSL_EINVAL, 0);
  42.     }
  43.  
  44.   if (mu != 0 && mu != 1)
  45.     {
  46.       GSL_ERROR_VAL ("mu must be 0 or 1", GSL_EINVAL, 0);
  47.     }
  48.  
  49.   if (nu != 0 && nu != 1)
  50.     {
  51.       GSL_ERROR_VAL ("nu must be 0 or 1", GSL_EINVAL, 0);
  52.     }
  53.  
  54.   t = (gsl_integration_qaws_table *) 
  55.     malloc(sizeof(gsl_integration_qaws_table));
  56.  
  57.   if (t == 0)
  58.     {
  59.       GSL_ERROR_VAL ("failed to allocate space for qaws_table struct",
  60.             GSL_ENOMEM, 0);
  61.     }
  62.  
  63.   t->alpha = alpha;
  64.   t->beta = beta;
  65.   t->mu = mu;
  66.   t->nu = nu;
  67.   
  68.   initialise (t->ri, t->rj, t->rg, t->rh, alpha, beta);
  69.   
  70.   return t;
  71. }
  72.  
  73.  
  74. int
  75. gsl_integration_qaws_table_set (gsl_integration_qaws_table * t,
  76.                                 double alpha, double beta, int mu, int nu)
  77. {
  78.   if (alpha < -1.0)
  79.     {
  80.       GSL_ERROR ("alpha must be greater than -1.0", GSL_EINVAL);
  81.     }
  82.  
  83.   if (beta < -1.0)
  84.     {
  85.       GSL_ERROR ("beta must be greater than -1.0", GSL_EINVAL);
  86.     }
  87.  
  88.   if (mu != 0 && mu != 1)
  89.     {
  90.       GSL_ERROR ("mu must be 0 or 1", GSL_EINVAL);
  91.     }
  92.  
  93.   if (nu != 0 && nu != 1)
  94.     {
  95.       GSL_ERROR ("nu must be 0 or 1", GSL_EINVAL);
  96.     }
  97.  
  98.   t->alpha = alpha;
  99.   t->beta = beta;
  100.   t->mu = mu;
  101.   t->nu = nu;
  102.   
  103.   initialise (t->ri, t->rj, t->rg, t->rh, alpha, beta);
  104.  
  105.   return GSL_SUCCESS;
  106. }
  107.  
  108.  
  109. void
  110. gsl_integration_qaws_table_free (gsl_integration_qaws_table * t)
  111. {
  112.   free (t);
  113. }
  114.  
  115. static void
  116. initialise (double * ri, double * rj, double * rg, double * rh,
  117.         double alpha, double beta)
  118. {
  119.   const double alpha_p1 = alpha + 1.0;
  120.   const double beta_p1 = beta + 1.0;
  121.  
  122.   const double alpha_p2 = alpha + 2.0;
  123.   const double beta_p2 = beta + 2.0;
  124.  
  125.   const double r_alpha = pow (2.0, alpha_p1);
  126.   const double r_beta = pow (2.0, beta_p1);
  127.  
  128.   size_t i;
  129.   
  130.   double an, anm1;
  131.  
  132.   ri[0] = r_alpha / alpha_p1;
  133.   ri[1] = ri[0] * alpha / alpha_p2;
  134.  
  135.   an = 2.0;
  136.   anm1 = 1.0;
  137.  
  138.   for (i = 2; i < 25; i++)
  139.     {
  140.       ri[i] = -(r_alpha + an * (an - alpha_p2) * ri[i - 1])
  141.     / (anm1 * (an + alpha_p1));
  142.       anm1 = an;
  143.       an = an + 1.0;
  144.     }
  145.  
  146.   rj[0] = r_beta / beta_p1;
  147.   rj[1] = rj[0] * beta / beta_p2;
  148.  
  149.   an = 2.0;
  150.   anm1 = 1.0;
  151.  
  152.   for (i = 2; i < 25; i++)
  153.     {
  154.       rj[i] = -(r_beta + an * (an - beta_p2) * rj[i - 1])
  155.     / (anm1 * (an + beta_p1));
  156.       anm1 = an;
  157.       an = an + 1.0;
  158.     }
  159.  
  160.   rg[0] = -ri[0] / alpha_p1;
  161.   rg[1] = -rg[0] - 2.0 * r_alpha / (alpha_p2 * alpha_p2);
  162.  
  163.   an = 2.0;
  164.   anm1 = 1.0;
  165.  
  166.   for (i = 2; i < 25; i++)
  167.     {
  168.       rg[i] = -(an * (an - alpha_p2) * rg[i - 1] - an * ri[i - 1]
  169.         + anm1 * ri[i]) / (anm1 * (an + alpha_p1));
  170.       anm1 = an;
  171.       an = an + 1.0;
  172.     }
  173.  
  174.   rh[0] = -rj[0] / beta_p1;
  175.   rh[1] = -rh[0] - 2.0 * r_beta / (beta_p2 * beta_p2);
  176.  
  177.   an = 2.0;
  178.   anm1 = 1.0;
  179.  
  180.   for (i = 2; i < 25; i++)
  181.     {
  182.       rh[i] = -(an * (an - beta_p2) * rh[i - 1] - an * rj[i - 1]
  183.         + anm1 * rj[i]) / (anm1 * (an + beta_p1));
  184.       anm1 = an;
  185.       an = an + 1.0;
  186.     }
  187.  
  188.   for (i = 1; i < 25; i += 2)
  189.     {
  190.       rj[i] *= -1;
  191.       rh[i] *= -1;
  192.     }
  193. }
  194.